# Wrangle the capture data
captures <- # extract the capture datatable captures <- # extract the capture
captures <- # extract the capture datatable datatable
dbReadTable(CeutaCLOSED, "Captures") %>%  dbReadTable(CeutaCLOSED, "Captures")
dbReadTable(CeutaCLOSED, "Captures") %>%  %>%
# subset to only include snowy plovers
dplyr::filter(species == "SNPL") %>%  dplyr::filter(species == "SNPL") %>%
# group by individual
group_by(ring) %>%  group_by(ring) %>%
# determine the first year that an individual was seen
mutate(originalsight = ifelse(age == "J", year, year[which.min(as.numeric(year))])) %>%
# determine if the first encouter was as an adult or a juvenile
mutate(chick = age[which.min(as.numeric(year))]) %>%  mutate(chick =
mutate(chick = age[which.min(as.numeric(year))]) %>%  age[which.min(as.numeric(year))])
mutate(chick = age[which.min(as.numeric(year))]) %>%  %>%
# specify first encounters as juveniles as recruits and adults as immigrants
mutate(recruit = ifelse(chick == "J", "Recruit", "Immigrant")) %>%  mutate(recruit
mutate(recruit = ifelse(chick == "J", "Recruit", "Immigrant")) %>%  =
mutate(recruit = ifelse(chick == "J", "Recruit", "Immigrant")) %>%  ifelse(chick
mutate(recruit = ifelse(chick == "J", "Recruit", "Immigrant")) %>%  == "J",
mutate(recruit = ifelse(chick == "J", "Recruit", "Immigrant")) %>%  "Recruit",
mutate(recruit = ifelse(chick == "J", "Recruit", "Immigrant")) %>%  "Immigrant"))
mutate(recruit = ifelse(chick == "J", "Recruit", "Immigrant")) %>%  %>%
# if a recruit then specify the year at which it was first encountered (i.e.,
# it's birth year); 0 for immigrants
mutate(birth = ifelse(recruit == "Recruit", originalsight, "0")) %>%  mutate(birth
mutate(birth = ifelse(recruit == "Recruit", originalsight, "0")) %>%  =
mutate(birth = ifelse(recruit == "Recruit", originalsight, "0")) %>%  ifelse(recruit
mutate(birth = ifelse(recruit == "Recruit", originalsight, "0")) %>%  ==
mutate(birth = ifelse(recruit == "Recruit", originalsight, "0")) %>%  "Recruit",
mutate(birth = ifelse(recruit == "Recruit", originalsight, "0")) %>%  originalsight,
mutate(birth = ifelse(recruit == "Recruit", originalsight, "0")) %>%  "0")) %>%
# remove 13 rows that do not have ring information (i.e., these are
# 'observation only' captures of XX.XX|XX.XX birds)
dplyr::filter(ring != "NA") %>%  dplyr::filter(ring != "NA") %>%
# subset to females
dplyr::filter(sex == "F")

# Assess the sample size of the capture population assess how many individuals
# we have: Female = 860
captures %>%
    group_by(sex) %>%
    summarise(count = n_distinct(ring)) %>%
    collect() %>%
    kable(col.names = c("Sex", "Number of individuals")) %>%
    kable_styling() %>%
    scroll_box(width = "50%")

# assess how many immigrants and recruits we have: Immigrant Females = 406,
# Recruit Females = 454
captures %>%
    group_by(recruit, sex) %>%
    summarise(n_inds = n_distinct(ring)) %>%
    collect() %>%
    kable(col.names = c("Status", "Sex", "Number of individuals")) %>%
    kable_styling() %>%
    scroll_box(width = "50%")

# Wrangle the resight data
resightings <- # extract the resight datatable resightings <- # extract the
resightings <- # extract the resight datatable resight datatable
dbReadTable(CeutaCLOSED, "Resights") %>%  dbReadTable(CeutaCLOSED, "Resights")
dbReadTable(CeutaCLOSED, "Resights") %>%  %>%
# subset to only include snowy plovers
dplyr::filter(species == "SNPL") %>%  dplyr::filter(species == "SNPL") %>%
# merge the ring identity to the resight codes
left_join(., dplyr::select(captures, ring, code), by = "code") %>%  left_join(.,
left_join(., dplyr::select(captures, ring, code), by = "code") %>%  dplyr::select(captures,
left_join(., dplyr::select(captures, ring, code), by = "code") %>%  ring,
left_join(., dplyr::select(captures, ring, code), by = "code") %>%  code), by =
left_join(., dplyr::select(captures, ring, code), by = "code") %>%  "code") %>%
# remove duplicates
distinct()

# extract combinations that are unique
unique_combos <- resightings %>%  unique_combos <- resightings %>%
# remove all combos without a distinct set of rings (i.e., needs at most 4 X's)
dplyr::filter(str_count(code, "X") < 5) %>%  dplyr::filter(str_count(code, "X")
dplyr::filter(str_count(code, "X") < 5) %>%  < 5) %>%
# remove all combos with uncertainty (i.e., no ?'s)
dplyr::filter(str_detect(code, "\?", negate = TRUE)) %>%  dplyr::filter(str_detect(code,
dplyr::filter(str_detect(code, "\?", negate = TRUE)) %>%  "\?", negate = TRUE))
dplyr::filter(str_detect(code, "\?", negate = TRUE)) %>%  %>%
# remove all combos without the appropriate number of rings
dplyr::filter(nchar(code) == 11) %>%  dplyr::filter(nchar(code) == 11) %>%
# extract the combos that have only one metal ring associate with them
group_by(code) %>%
    summarise(n_obs = n_distinct(ring)) %>%
    arrange(desc(n_obs)) %>%
    dplyr::filter(n_obs == 1)

# Tidy up capture and resight data
resightings <- resightings %>%  resightings <- resightings %>%
# subset resightings to only include observations of unique combinations
dplyr::filter(code %in% unique_combos$code) %>%  dplyr::filter(code %in%
dplyr::filter(code %in% unique_combos$code) %>%  unique_combos$code) %>%
# remove all combos that don't have a ring associated with them
dplyr::filter(!is.na(ring)) %>%  dplyr::filter(!is.na(ring)) %>%
# sort by ring and year to assess the output
arrange(ring, year)

capture_final <- captures %>%  capture_final <- captures %>%
# paste ring and year together to make a unique identifier for this observation
unite(ring_year, ring, year, remove = FALSE) %>%  unite(ring_year, ring, year,
unite(ring_year, ring, year, remove = FALSE) %>%  remove = FALSE) %>%
# specify as a capture
mutate(observation = "capture") %>%  mutate(observation = "capture") %>%
# clean up output
dplyr::select(ring_year, sex, ring, year, recruit, birth, observation) %>%
    arrange(ring, year)

resightings_final <- resightings %>%  resightings_final <- resightings %>%
# paste ring and year together to make a unique identifier for this observation
unite(ring_year, ring, year, remove = FALSE) %>%  unite(ring_year, ring, year,
unite(ring_year, ring, year, remove = FALSE) %>%  remove = FALSE) %>%
# specify as a resight
mutate(observation = "resight") %>%  mutate(observation = "resight") %>%
# clean up output
dplyr::select(ring_year, ring, year, observation) %>%  dplyr::select(ring_year,
dplyr::select(ring_year, ring, year, observation) %>%  ring, year, observation)
dplyr::select(ring_year, ring, year, observation) %>%  %>%
# join with the capture data to merge recruit status, sex, and birth year
left_join(., dplyr::select(capture_final, ring, recruit, birth, sex), by = "ring") %>%
# remove duplicates
distinct()

# Bind capture and resight data into a complete encounter history
encounter_histories <- bind_rows(capture_final, resightings_final) %>%
    arrange(ring_year)

# Remove resight encounters that occurred prior to the first capture determine
# the year of first capture for each individual
first_cap <- capture_final %>%
    group_by(ring) %>%
    dplyr::filter(as.numeric(year) == min(as.numeric(year))) %>%
    dplyr::select(ring, year) %>%
    distinct(ring, .keep_all = TRUE) %>%
    rename(first_cap = year) %>%
    arrange(first_cap)

# determine which resights occurred before the first capture
resights_before_first_capture <- encounter_histories %>%
    left_join(., first_cap, by = "ring") %>%
    dplyr::filter(observation == "resight" & (as.numeric(year) < as.numeric(first_cap)))

# exclude early resightings from encounter history
encounter_histories <- encounter_histories %>%
    dplyr::filter(ring_year %!in% resights_before_first_capture$ring_year) %>%
    arrange(ring, as.numeric(year))

# Make the encounter history table needed for the survival analysis
encounter_history_table <- distinct(encounter_histories, ring_year, .keep_all = TRUE) %>%
    dplyr::select(ring, year) %>%
    arrange(ring) %>%
    mutate(year = as.integer(year)) %>%
    CensusToCaptHist(ID = .$ring, d = .$year, timeInt = "Y") %>%
    mutate(ring = rownames(.), ID = as.character(ID))

# Extract the known birth and death information for each individual
birth_death_mat <- encounter_histories %>%
    dplyr::select(ring, birth) %>%
    mutate(death = 0) %>%
    arrange(ring) %>%
    distinct(ring, .keep_all = TRUE)

# Unite the encounter history table with the birth and death matrices
raw_life_table_females_2006_2020 <- left_join(birth_death_mat, encounter_history_table,
    by = "ring") %>%
    as.data.frame() %>%
    distinct() %>%
    dplyr::select(ID, birth, death, "2006", "2007", "2008", "2009", "2010", "2011",
        "2012", "2013", "2014", "2015", "2016", "2017", "2018", "2019", "2020", ring) %>%
    mutate_at(vars(-ring), as.numeric) %>%
    mutate(sum_years = rowSums(.[4:18])) %>%
    dplyr::filter(sum_years > 1 | birth == 0) %>%
    mutate(ID = as.numeric(row.names(.))) %>%
    dplyr::select(-sum_years)

# Run 'DataCheck' BaSTA function to make final cleans to encounter history. The
# only change needed is that the birth year of recruits should be '0' instead
# of '1'; this function will solve this issue and provide a cleaned version
# ready for analysis. In total, we have 918 encounters of 452 individually
# marked females that were seen as adults, of which 46 are of known birth year.
BaSTA_checked_life_table_females_2006_2020 <- DataCheck(object = raw_life_table_females_2006_2020[,
    -c(length(raw_life_table_females_2006_2020))], studyStart = 2006, studyEnd = 2020,
    autofix = rep(1, 7), silent = FALSE)

load(file = "data/BaSTA_checked_life_table_females_2006-2020.rds")

multiout_females <- multibasta(object = BaSTA_checked_life_table_females_2006_2020$newData,
    studyStart = 2006, studyEnd = 2020, covarsStruct = "fused", minAge = 1, niter = 8e+05,
    burnin = 10000, thinning = 2000, nsim = 4, parallel = TRUE, ncpus = 4, updateJumps = TRUE)

load(file = "output/multibasta_output_females_min_age_1_2006-2020.rds")
load(file = "data/raw_life_table_females_2006_2020.rds")

# DIC model selection shows that the Logistic mortality model with bathtub 
# shape fits our data the best
BaSTA_DIC_table <- 
  multiout_females$DICs %>% 
  as.data.frame() %>% 
  mutate(model = ifelse(model == "LO", "Logistic",
                        ifelse(model == "WE", "Weibull",
                               ifelse(model == "GO", "Gompertz", "Exponential"))),
         shape = paste0(toupper(substr(shape, 1, 1)), 
                        substr(shape, 2, nchar(.)))) %>% 
  mutate(shape = ifelse(shape == "Simp", "Simple", shape)) %>% 
  dplyr::select(model, shape, k, DICdiff) %>% 
  gt() %>% 
  cols_label(model = "Mortality function",
             shape = "Shape",
             k = md("*k*"),
             # DIC = md("*DIC*"),
             DICdiff = md("\U0394*DIC*")) %>% 
  fmt_number(columns = vars(DICdiff),
             decimals = 2,
             use_seps = FALSE) %>% 
  # tab_options(column_labels.font.weight = "bold",
  #             table.width = pct(50),
  #             column_labels.font.size = 14,
  #             table.font.size = 12,
  #             data_row.padding = 5, table.align = "left") %>% 
  fmt_missing(columns = 3,
              missing_text = "") %>% 
  tab_options(row_group.font.weight = "bold",
              row_group.background.color = brewer.pal(9,"Greys")[3],
              table.font.size = 12,
              data_row.padding = 3,
              row_group.padding = 4,
              summary_row.padding = 2,
              column_labels.font.size = 14,
              row_group.font.size = 12,
              table.width = pct(60))

BaSTA_DIC_table

load("output/BaSTA_age_estimates_2006-2020.rds")

#### extract females and their nests ####
nest_caps_F <-
  
  # extract snowy plover captures of adult females
  dbReadTable(CeutaCLOSED,"Captures") %>% 
  dplyr::filter(species == "SNPL" & sex == "F" & age == "A") %>% 
  
  # summarise body morphometric data with the earliest capture of each nest
  group_by(ID, ring, code, sex, age) %>% 
  summarise(ad_cap_date = date[which.min(date)],
            ad_cap_time = time[which.min(date)],
            ad_weight = weight[which.min(date)],
            left_tarsus = mean(left_tarsus, na.rm = TRUE),
            right_tarsus = mean(right_tarsus, na.rm = TRUE),
            ad_bill = mean(bill, na.rm = TRUE),
            ad_fat = fat[which.min(date)]) %>% 
  mutate(ad_tarsi = mean(c(left_tarsus, right_tarsus), na.rm = TRUE)) %>% 
  group_by(ring) %>% 
  mutate(avg_ad_tarsi = mean(ad_tarsi, na.rm = TRUE),
         sd_ad_tarsi = sd(ad_tarsi, na.rm = TRUE),
         n_ad_tarsi = sum(!is.na(ad_tarsi)),
         avg_ad_weight = mean(ad_weight, na.rm = TRUE),
         sd_ad_weight = sd(ad_weight, na.rm = TRUE),
         n_ad_weight = sum(!is.na(ad_weight))) %>% 

  # join the nest data
  left_join(., dplyr::filter(dbReadTable(CeutaCLOSED, "Nests"), species == "SNPL"), 
            by = "ID") %>% 
  
  # remove repeated observations of nests 2020_C_4 and 2020_D_201 that were done
  # for training purposes in the field
  dplyr::filter(ID != "2020_C_4" | observer != "DVR") %>% 
  dplyr::filter(ID != "2020_D_201" | observer != "MC") %>%
  
  # remove duplicated rows
  distinct() %>% 
  
  # set dates as numeric
  mutate(nest_initiation_date = as.numeric(nest_initiation_date),
         end_date = as.numeric(end_date),
         found_date = as.numeric(found_date),
         year = as.factor(year)) %>% 
  
  # change to as.Date format
  plover_date_convert(input = "Rdate") %>%
  
  group_by(ID) %>% 
  # specify the lay date as the nest initiation date. If this is unknown, then
  # subtract 25 days from the end date if the nest hatched. If this is not the case,
  # then subtract 11 days from the found date if the float score of the first 
  # egg is "F". If this isn't the case, then take the found date.
  mutate(lay_date1 = 
           as.Date(
                      ifelse((!is.na(found_date) & float1 == "A"), found_date - 0,
                        ifelse((!is.na(found_date) & float1 == "AB"), found_date - 1,
                          ifelse((!is.na(found_date) & float1 == "B"), found_date - 2,
                            ifelse((!is.na(found_date) & float1 == "C"), found_date - 5,
                              ifelse((!is.na(found_date) & float1 == "D"), found_date - 8,
                                ifelse((!is.na(found_date) & float1 == "E"), found_date - 10,
                                  ifelse((!is.na(found_date) & float1 == "F"), found_date - 11,
                                    ifelse((!is.na(end_date) & fate == "Hatch"), end_date - 25, 
                                      ifelse(!is.na(found_date), found_date - 17, NA))))))))), 
             origin = "1970-01-01"),
         lay_date2 = 
           as.Date(
                 ifelse((!is.na(found_date) & float2 == "A"), found_date - 0,
                   ifelse((!is.na(found_date) & float2 == "AB"), found_date - 1,
                     ifelse((!is.na(found_date) & float2 == "B"), found_date - 2,
                       ifelse((!is.na(found_date) & float2 == "C"), found_date - 5,
                         ifelse((!is.na(found_date) & float2 == "D"), found_date - 8,
                           ifelse((!is.na(found_date) & float2 == "E"), found_date - 10,
                             ifelse((!is.na(found_date) & float2 == "F"), found_date - 11,
                                ifelse((!is.na(end_date) & fate == "Hatch"), end_date - 25, 
                                  ifelse(!is.na(found_date), found_date - 17, NA))))))))), 
             origin = "1970-01-01"),
         lay_date3 = 
           as.Date(
                 ifelse((!is.na(found_date) & float3 == "A"), found_date - 0,
                   ifelse((!is.na(found_date) & float3 == "AB"), found_date - 1,
                     ifelse((!is.na(found_date) & float3 == "B"), found_date - 2,
                       ifelse((!is.na(found_date) & float3 == "C"), found_date - 5,
                         ifelse((!is.na(found_date) & float3 == "D"), found_date - 8,
                           ifelse((!is.na(found_date) & float3 == "E"), found_date - 10,
                              ifelse((!is.na(found_date) & float3 == "F"), found_date - 11,
                                 ifelse((!is.na(end_date) & fate == "Hatch"), end_date - 25, 
                                   ifelse(!is.na(found_date), found_date - 17, NA))))))))), 
             origin = "1970-01-01"),
         min_egg_float = min(c(float1, float2, float3), na.rm = TRUE)) %>% 
  
  mutate(lay_date = as.Date(min(c(lay_date1, lay_date3, lay_date3), 
                                na.rm = TRUE), 
                            origin = "1970-01-01")) %>% 
  mutate(lay_date = as.Date(ifelse(min_egg_float == "F" & !is.na(found_date) & fate != "Hatch", 
                                   found_date - 17, 
                                   ifelse((min_egg_float == "F" | is.na(min_egg_float)) & !is.na(found_date) & fate == "Hatch", 
                                          end_date - 25, lay_date)), origin = "1970-01-01")) %>% 
  mutate(lay_date = as.Date(ifelse(is.na(lay_date) & !is.na(found_date), 
                                   found_date - 17, lay_date), 
                            origin = "1970-01-01")) %>% 

  mutate(lay_date_method =
           ifelse((min_egg_float < "F" & !is.na(min_egg_float)), "floatation",
                    ifelse((!is.na(end_date) & fate == "Hatch"), "hatch date",
                                  ifelse(!is.na(found_date), "found_date - 17", NA)))) %>%
  
  # select(ID, fate, found_date, end_date, lay_date1, lay_date2, lay_date3,
  #        float1, float2, float3, min_egg_float, lay_date, lay_date_method) %>%
  # View()
  
  # create a julian lay date
  mutate(jul_lay_date = as.numeric(format(lay_date, "%j"))) %>% 
  
  # remove any duplicated rows resulting from joining above
  distinct() %>% 
  
  # mutate(
  #   # scale the julian lay date by year
  #   jul_std_date = scale_by(jul_lay_date ~ year, ., scale = 0)) %>%
  
  # remove any rows without lay date information
  dplyr::filter(!is.na(jul_lay_date)) %>% 
  
  # remove rows with NA egg dimensions
  dplyr::filter(!is.na(width1) & !is.na(length1))

# determine how many years of observations each female has
n_years <- 
  nest_caps_F %>% 
  group_by(ring) %>% 
  summarise(n_years_obs = n_distinct(year))

# join repeated meausres info and consolidate dataframe
nest_caps_F <-
  left_join(nest_caps_F, n_years, by = "ring") %>% 
  dplyr::select(-left_tarsus, -right_tarsus, -species, -population, 
                -site, -nest, -found_date, -found_time, -nest_initiation_date, 
                -photo, -observer, -comments)

# extract dimensions of egg 1
nest_caps_F_egg1 <- 
  nest_caps_F %>% 
  dplyr::select(-width2, -width3, 
                -length2, -length3,
                -float2, -float3) %>%
  rename(width = width1,
         length = length1,
         float = float1)

# extract dimensions of egg 2
nest_caps_F_egg2 <- 
  nest_caps_F %>% 
  dplyr::select(-width1, -width3, 
                -length1, -length3,
                -float1, -float3) %>%
  rename(width = width2,
         length = length2,
         float = float2)

# extract dimensions of egg 3
nest_caps_F_egg3 <- 
  nest_caps_F %>% 
  dplyr::select(-width2, -width1, 
                -length2, -length1,
                -float2, -float1) %>%
  rename(width = width3,
         length = length3,
         float = float3)

# bind all egg measurements and remove NA observations
nest_caps_F <- 
  bind_rows(nest_caps_F_egg1, 
            nest_caps_F_egg2, 
            nest_caps_F_egg3) %>% 
  dplyr::filter(!is.na(width) | !is.na(length)) %>% 
  
  # calculate egg volume
  mutate(volume = 0.486 * length * width^2,
         year = as.integer(as.character(year))) %>% 
  
  # transform measures to cm
  mutate(length_cm = length/10,
         width_cm = width/10,
         volume_cm = volume/1000) %>% 
  
  # sort by ring and nest ID and remove duplicated rows
  arrange(ring, ID) %>% 
  distinct() %>% 
  
  # merge with basta age estimates
  left_join(., BaSTA_ages, by = "ring") %>% 
  
  mutate(
    # calculate the estimated age at a given year
    est_age = year - est_b,
    
    # calculate upper and lower 95% CI for age at a given year
    est_age_lower = year - upper_b,
    est_age_upper = year - lower_b,
    
    est_death_age = upper_d) %>% 
  
  # remove duplicate rows and extraneous columns
  distinct()

# extract egg laying order from float scores
first_egg <- 
  nest_caps_F %>% 
  arrange(ID, desc(float)) %>% 
  group_by(ID) %>% 
  slice(1) %>% 
  mutate(egg = "egg1")

second_egg <-
  nest_caps_F %>% 
  arrange(ID, desc(float)) %>% 
  group_by(ID) %>% 
  slice(2) %>%
  mutate(egg = "egg2")

third_egg <-
  nest_caps_F %>% 
  arrange(ID, desc(float)) %>% 
  group_by(ID) %>% 
  slice(3) %>%
  mutate(egg = "egg3")

nest_caps_F <- 
  bind_rows(first_egg, 
            second_egg, 
            third_egg)

distinct_first_eggs <- 
  nest_caps_F %>% 
  dplyr::filter(!is.na(float)) %>% 
  dplyr::filter(egg %in% c("egg1", "egg2")) %>% 
  group_by(ID) %>% 
  summarise(n_float_scores = n_distinct(float),
            first_egg = nth(float, 1),
            second_egg = nth(float, 2)) %>% 
  arrange(desc(n_float_scores))

nest_caps_F <- 
  left_join(nest_caps_F, dplyr::select(distinct_first_eggs, 
                                       ID, n_float_scores), by = "ID") %>% 
  rename(old_egg1 = n_float_scores) %>% 
  mutate(old_egg1 = ifelse(old_egg1 == 2, 1, 0))

# extract the first, second, third, and fourth nests for each individual in a
# given year
first_nests <- 
  nest_caps_F %>% 
  dplyr::select(ring, year, ID, jul_lay_date) %>%
  arrange(ring, year, jul_lay_date) %>%
  distinct() %>% 
  group_by(ring, year) %>% 
  slice(1) %>% 
  mutate(nest_order = 1) %>% 
  dplyr::select(-jul_lay_date) %>% 
  mutate(ring_year = paste(ring, year, sep = "_"))

second_nests <- 
  nest_caps_F %>% 
  dplyr::select(ring, year, ID, jul_lay_date) %>%
  arrange(ring, year, jul_lay_date) %>%
  distinct() %>% 
  group_by(ring, year) %>% 
  slice(2) %>% 
  mutate(nest_order = 2) %>% 
  dplyr::select(-jul_lay_date) %>% 
  mutate(ring_year = paste(ring, year, sep = "_"))

third_nests <- 
  nest_caps_F %>% 
  dplyr::select(ring, year, ID, jul_lay_date) %>%
  arrange(ring, year, jul_lay_date) %>%
  distinct() %>% 
  group_by(ring, year) %>% 
  slice(3) %>% 
  mutate(nest_order = 3) %>% 
  dplyr::select(-jul_lay_date) %>% 
  mutate(ring_year = paste(ring, year, sep = "_"))

fourth_nests <- 
  nest_caps_F %>% 
  dplyr::select(ring, year, ID, jul_lay_date) %>%
  arrange(ring, year, jul_lay_date) %>%
  distinct() %>% 
  group_by(ring, year) %>% 
  slice(4) %>% 
  mutate(nest_order = 4) %>% 
  dplyr::select(-jul_lay_date) %>% 
  mutate(ring_year = paste(ring, year, sep = "_"))

# bind all together into a single dataframe
nest_caps_F <- 
  bind_rows(first_nests, second_nests, third_nests, fourth_nests) %>% 
  left_join(nest_caps_F, ., by = c("ring", "year", "ID"))

# determine the fate of the first nest of the season for each female
nest_caps_F <- 
  nest_caps_F %>% 
  group_by(ring_year) %>% 
  summarise(nest_1_fate = fate[which.min(jul_lay_date)]) %>% 
  left_join(nest_caps_F, ., by = "ring_year")

# Determine the age at first capture for females in the data
age_at_first_cap_info <- 
  dbReadTable(CeutaCLOSED, "Captures") %>%
  dplyr::filter(ring %in% nest_caps_F$ring) %>%
  plover_date_convert(input = "Rdate") %>%
  group_by(ring) %>%
  summarise(age_first_cap = age[which.min(year)],
            year_first_cap = as.numeric(year[which.min(year)])) %>%
  collect() %>% 
  
  # manually assign J to CA2036 and CA1526 which were captured as chicks during 
  # the small study in 2004
  mutate(age_first_cap = ifelse(ring %in% c("CA2036", "CA1526"), "J", age_first_cap))

# join the age at first capture info to the data
nest_caps_F <- 
  left_join(nest_caps_F, age_at_first_cap_info, 
            by = "ring") %>% 
  mutate(year_first_cap = ifelse(age_first_cap == "J", est_b, year_first_cap)) %>% 
  mutate(conservative_age = year - year_first_cap)

# Need to check these four nests
nest_caps_F %>% 
  dplyr::filter(code != female) %>%
  dplyr::select(year, ID, male, female, ring, code) %>%  
  distinct
nest_caps_F %>% 
  dplyr::filter(code == male) %>%
  dplyr::select(year, ID, male, female, ring, code) %>%  
  distinct

# check tarsus measurements
nest_caps_F %>% 
  # dplyr::filter(sd_ad_tarsi > 1) %>% 
  dplyr::select(ring, avg_ad_tarsi, n_ad_tarsi, sd_ad_tarsi) %>% 
  distinct() %>% 
  arrange(desc(sd_ad_tarsi))

# extract polyandry and multiclutching information
nest_caps_F <-
  nest_caps_F %>% 
  dplyr::select(year, ID, male, female, ring) %>%
  arrange(female) %>%
  group_by(female, ring, year) %>%
  summarise(n_mates = n_distinct(male, na.rm = TRUE),
            n_nests = n_distinct(ID, na.rm = TRUE)) %>% 
  mutate(polyandry = ifelse(n_mates > 1, "poly", "mono"),
         multiclutch = ifelse(n_nests > 1, "multi", "single"),
         year = as.integer(year),
         n_mates = ifelse(n_mates == 0, 1, n_mates)) %>% 
  left_join(nest_caps_F, ., 
            by = c("year", "ring", "female")) %>% 
  distinct

# add average, first, and last age information to each ring
age_summary <- 
  function(df){
    
    # extract the average, first, and last age for each individual
    ring_Age <- 
      df %>%
      dplyr::select(ring, est_age, conservative_age, 
                    est_age_lower, est_age_upper) %>%
      distinct() %>% 
      group_by(ring) %>% 
      summarise(first_age_t = min(est_age) - 1,
                conservative_first_age = min(conservative_age),
                first_age_lower = min(est_age_lower),
                first_age_upper = min(est_age_upper),
                last_age_t = max(est_age) - 1,
                conservative_last_age = max(conservative_age),
                last_age_lower = max(est_age_lower),
                last_age_upper = max(est_age_upper))
    
    # merge with dataframe
    df2 <- 
      left_join(df, ring_Age, by = "ring") %>% 
      mutate(est_age_t = est_age - 1,
             conservative_age = ifelse(age_first_cap == "J", 
                                       conservative_age - 1, 
                                       conservative_age),
             conservative_first_age_t = ifelse(age_first_cap == "J", 
                                               conservative_first_age - 1, 
                                               conservative_first_age),
             conservative_last_age_t = ifelse(age_first_cap == "J", 
                                              conservative_last_age - 1, 
                                              conservative_last_age)) %>% 
      mutate(est_age_t_deviation = est_age_t - first_age_t)
    
    
    return(df2)
    
  }

# apply age summary to data
nest_caps_F <-
  nest_caps_F %>% 
  age_summary(df = .) %>% 
  mutate(year = as.factor(year),
         ID = as.factor(ID),
         ring = as.factor(ring)) %>% 
  distinct()

eggs_2006_2020 <-
  nest_caps_F %>% 
  ungroup() %>% 
  mutate(ID = as.character(ID),
         ring = as.character(ring)) %>%
  
  # remove CN0424 since sex is unclear
  dplyr::filter(ring != "CN0424") %>% 
  
  # remove outlier egg measurements below 6000 cubic mm
  dplyr::filter(volume > 6000) #%>%

eggs_2006_2020 <- 
  eggs_2006_2020 %>% 
  dplyr::select(ID, jul_lay_date, year, ring_year) %>% 
  distinct() %>% 
  mutate(
    # scale the julian lay date by year
    jul_lay_date_std = scale_by(jul_lay_date ~ year, ., scale = 0)) %>% 
  
  # make the scaled date variable numeric class
  mutate(jul_lay_date_std_num = as.numeric(jul_lay_date_std)) %>% 
  group_by(ring_year) %>% 
  arrange(ring_year, jul_lay_date) %>% 
  mutate(laydate_deviation = jul_lay_date_std_num - jul_lay_date_std_num[which.min(jul_lay_date_std_num)],
         first_laydate = jul_lay_date_std_num[which.min(jul_lay_date_std_num)],
         last_laydate = jul_lay_date_std_num[which.max(jul_lay_date_std_num)]) %>% 
  left_join(eggs_2006_2020, ., 
            by = c("ID", "jul_lay_date", "year", "ring_year")) %>% 
  distinct()

#### Chick data wrangle ----
# Extract chick measurements from the capture data
chicks_measurements <-
  # read the Captures table
  dbReadTable(CeutaCLOSED, "Captures") %>% 
  # subset to SNPL
  dplyr::filter(species == "SNPL") %>% 
  # average the left and right tarsus length measurements in to one value
  mutate(tarsus = rowMeans(cbind(as.numeric(left_tarsus), 
                                 as.numeric(right_tarsus)), na.rm = TRUE)) %>% 
  # calculate Body Mass Index based on weight and structural size
  mutate(BMI = weight/(tarsus^2)) %>% 
  # subset to juvenile captures
  dplyr::filter(age == "J",
                # remove observations that are missing information for weight and tarsus
                !is.na(weight) & !is.na(tarsus) & !is.na(date)) %>% 
  # convert date columns to the %Y-%m-%d format
  plover_date_convert(input = "Rdate") %>% 
  # group by bird identity
  group_by(ring) %>%
  # specify the biological nest ID as the ID of the earliest capture (i.e.,
  # brood mixing can occur and create multiple nest IDs for a chick in the
  # capture data)
  dplyr::mutate(bio_ID = ID[which.min(date)]) %>% 
  # filter(bio_ID != ID) %>% 
  # convert to dataframe
  data.frame() %>% 
  # select the relevant columns
  dplyr::select(ring, year, bio_ID, age, sex, date, time, weight, tarsus, bill, BMI) %>% 
  # extract only chicks that are in the egg data
  dplyr::filter(bio_ID %in% eggs_2006_2020$ID) %>% 
  arrange(ring)

# Extract hatch dates
hatch_dates <-
  # read the Nests table
  dbReadTable(CeutaCLOSED, "Nests") %>%
  # extract only the nests that contain chicks in the previous capture subset
  dplyr::filter(ID %in% chicks_measurements$bio_ID) %>%
  # subset the nests that have hatch date information
  dplyr::filter(fate == "Hatch") %>%
  # define as a dataframe
  data.frame() %>%
  # classify date columns in the appropriate format
  plover_date_convert(input = "Rdate") %>% 
  # subset the result as simply the nest ID and their respective hatch dates
  dplyr::select(ID, end_date) %>%
  # rename the columns
  rename(bio_ID = ID,
         hatch_date = end_date)

# combine chick measurements with hatch dates
chicks_2006_2020 <- 
  # join the hatch dates to the chick captures
  left_join(x = chicks_measurements, y = hatch_dates, by = "bio_ID") %>%
  # classify variables as factor or numeric and calculate age at capture
  ungroup() %>% 
  mutate(age = as.numeric(date - hatch_date),
         year = as.factor(year),
         weight = as.numeric(weight),
         tarsus = as.numeric(tarsus),
         ring = as.factor(ring),
         sex = as.factor(sex)) %>% 
  rename(chick_ring = ring,
         ID = bio_ID) %>% 
  # specify ages less than 0 as 0 (i.e., hatch dates represent the average hatch
  # date of a brood and thus the earliest chick to hatch in a nest could be up
  # 2 days earlier than the nest's hatch date)
  mutate(age = ifelse(age < 0, 0, age)) %>% 
  # scale the numeric variables in preparation for modeling
  mutate(age.z = scale(age),
         weight.z = scale(weight),
         tarsus.z = scale(tarsus),
         BMI.z = scale(BMI)) %>% 
  # remove individuals that don't have an age value (i.e., hatch date was NA)
  dplyr::filter(!is.na(age)) %>% 
  left_join(., dplyr::select(eggs_2006_2020, ID, ring)) %>% 
  rename(mother_ring = ring)

# summarize chick morphometric data by nest and filter to only include chicks
# measured at age 0 or 1
chick_size_summary <- 
  chicks_2006_2020 %>% 
  arrange(chick_ring, age) %>% 
  group_by(chick_ring) %>% 
  slice(1) %>% 
  dplyr::filter(age %in% c(0, 1)) %>% 
  group_by(ID) %>% 
  summarise(avg_chick_tarsus = mean(tarsus, na.rm = TRUE),
            sd_chick_tarsus = sd(tarsus, na.rm = TRUE),
            avg_chick_bill = mean(bill, na.rm = TRUE),
            sd_chick_bill = sd(bill, na.rm = TRUE),
            avg_chick_weight = mean(weight, na.rm = TRUE),
            sd_chick_weight = sd(weight, na.rm = TRUE),
            avg_chick_BMI = mean(BMI, na.rm = TRUE),
            sd_chick_BMI = sd(BMI, na.rm = TRUE)) %>% 
  
  # remove questionable outlier chick measurements (hatch date likely incorrect)
  dplyr::filter(ID != "2020_H_6")

# join the summarised data together by nest
ceuta_egg_chick_female_data <- 
  left_join(eggs_2006_2020, chick_size_summary, by = "ID")

# saveRDS(ceuta_egg_chick_female_data,
#      file = "data/Ceuta_egg_chick_female_data.rds")

# extract the first and last observations of each individual
ceuta_egg_chick_female_data_3_year <-  
  ceuta_egg_chick_female_data %>% 
  dplyr::filter(n_years_obs > 2)

max_min_obs <-
  encounter_histories %>% 
  dplyr::filter(ring %in% ceuta_egg_chick_female_data_3_year$ring) %>% 
  left_join(., dplyr::select(BaSTA_ages, ring, est_b), by = "ring") %>%
  mutate(est_age = as.numeric(year) - as.numeric(est_b)) %>% 
  group_by(ring) %>% 
  summarise(min_age = min(est_age),
            max_age = max(est_age),
            obs_span = max(est_age) - min(est_age)) %>% 
  as.data.frame() %>% 
  left_join(., dplyr::select(ceuta_egg_chick_female_data_3_year, ring, age_first_cap)) %>% 
  distinct()

# manually adjust the ages of CA1526 and CA2036 to reflect their first encounter
# during the 2004 pilot study as chicks (i.e., aged 0)
max_min_obs[which(max_min_obs$ring %in% c("CA1526", "CA2036")), "min_age"] <-
  c(0, 0)

# manually adjust the min age estimate of CA1579 to relect that it was first
# encountered as an adult during the 2004 pilot study (i.e., subtract 2 years
# from BaSTA estimate of age at first encounter)
max_min_obs[which(max_min_obs$ring == "CA1579"), "min_age"] <-
  max_min_obs[which(max_min_obs$ring == "CA1579"), "min_age"] - 2
max_min_obs[which(max_min_obs$ring == "CA1579"), "obs_span"] <-
  max_min_obs[which(max_min_obs$ring == "CA1579"), "max_age"] -
  max_min_obs[which(max_min_obs$ring == "CA1579"), "min_age"]

# specify the factor levels according to the rank (for plotting)
max_min_obs$ring_ordered <- 
  factor(max_min_obs$ring, 
         levels = unique(max_min_obs$ring[order(max_min_obs$min_age, 
                                                max_min_obs$obs_span, 
                                                decreasing = TRUE)]), 
         ordered = TRUE)

# calculate the number of clutches per year for each individual and join back to 
# ranked dataframe
eggdf_clutches <- 
  ceuta_egg_chick_female_data_3_year %>% 
  dplyr::group_by(ring, year) %>%
  dplyr::summarise(n_clutches = n_distinct(ID)) %>% 
  dplyr::right_join(dplyr::select(max_min_obs, 
                                  c(ring, min_age, obs_span, age_first_cap)), 
                    by = "ring") %>% 
  dplyr::right_join(dplyr::select(ceuta_egg_chick_female_data_3_year, c(ring, year, est_age, polyandry)), 
                    by = c("ring", "year")) %>% 
  distinct() %>% 
  mutate(polyandry = ifelse(polyandry == "poly", "2", "1"))

# specify the factor levels according to the rank (for plotting)
eggdf_clutches$ring_ordered <- 
  factor(eggdf_clutches$ring, 
         levels = unique(eggdf_clutches$ring[order(eggdf_clutches$min_age, 
                                                   eggdf_clutches$obs_span, 
                                                   decreasing = TRUE)]), 
         ordered = TRUE)

# create a ranking variable and confidence interval limits of age estimate
eggdf_n_upper <-
  ceuta_egg_chick_female_data_3_year %>%
  dplyr::group_by(ring) %>%
  dplyr::summarise(CI_age = max(est_death_age))


eggdf_age_CI <-
  ceuta_egg_chick_female_data_3_year %>%
  dplyr::group_by(ring) %>%
  dplyr::summarise(CI_age = 1) %>%
  dplyr::bind_rows(eggdf_n_upper) %>%
  mutate(ring = as.factor(ring)) %>% 
  arrange(ring) %>% 
  dplyr::left_join(., dplyr::select(max_min_obs, 
                                    c(ring, min_age, obs_span, age_first_cap)))

# specify the factor levels according to the rank (for plotting)
eggdf_age_CI$ring_ordered <- 
  factor(eggdf_age_CI$ring, 
         levels = unique(eggdf_age_CI$ring[order(eggdf_age_CI$min_age, 
                                                 eggdf_age_CI$obs_span,
                                                 decreasing = TRUE)]), 
         ordered = TRUE)

# # join polyandry information
# eggdf_mating_sys_all <-
#   dbReadTable(Ceuta_OPEN, "BirdRef") %>% 
#   dplyr::select(year, ID, male, female) %>%
#   filter(female %in% eggdf$ring) %>% 
#   arrange(female) %>%
#   group_by(female, year) %>%
#   summarise(n_mates = n_distinct(male, na.rm = TRUE)) %>% 
#   mutate(polyandry = ifelse(n_mates > 1, "poly", "mono"),
#          year = as.integer(year),
#          n_mates = ifelse(n_mates == 0, 1, n_mates)) %>% 
#   rename(ring = female)
# 
# matings_without_egg_obs <- 
#   anti_join(eggdf_mating_sys_all, eggdf, 
#             by = c("ring", "year")) %>% 
#   left_join(., dplyr::select(eggdf, ring, age_first_cap, est_b), by = "ring") %>% 
#   mutate(est_age = year - est_b) %>% 
#   distinct()

# # specify the factor levels according to the rank (for plotting)
# matings_without_egg_obs_plot <- 
#   matings_without_egg_obs %>% 
#   dplyr::filter(ring %in% ceuta_egg_chick_female_data_3_year$ring)
# 
# matings_without_egg_obs_plot$ring_ordered <- 
#   factor(matings_without_egg_obs_plot$ring, 
#          levels = levels(eggdf_age_CI$ring_ordered),
#          ordered = TRUE)
# 
# dplyr::filter(ceuta_egg_chick_female_data_3_year, is.na(volume) & !is.na(polyandry))

# specify point size for number of clutches per year
point_size <- c(1, 2, 3, 4)
point_colors <- c("black", "#f03b20")

# create the first captured as adult plot
Imm_plot <-
  ggplot2::ggplot() +
  geom_line(data = dplyr::filter(eggdf_age_CI, age_first_cap == "A"),
            aes(x = CI_age, y = ring_ordered), color = "grey90",
            size = 2.5, lineend = "round") +
  geom_linerange(data = dplyr::filter(max_min_obs, age_first_cap == "A"),
                 orientation = "y", aes(y = ring_ordered, x = min_age,
                                        xmin = min_age, xmax = max_age),
                 size = 0.85, color = "grey60") +
  geom_point(data = dplyr::filter(eggdf_clutches, age_first_cap == "A"), 
             aes(x = est_age, y = ring_ordered, size = as.factor(n_clutches), 
                 fill = as.factor(polyandry)), shape = 21) +
  # geom_point(data = dplyr::filter(matings_without_egg_obs_plot, age_first_cap == "A"),
  #            aes(x = est_age, y = ring_ordered), size = 2, shape = 4) +
  scale_x_continuous(limits = c(0, 24),
                     breaks = seq(0, 24, by = 2)) +
  ylab("First captured as adults") +
  xlab("Age ± 95% CrI") +
  luke_theme +
  theme(legend.position = c(0.75, 0.8),
        legend.title = element_text(size = 9),
        legend.text = element_text(size = 9)) +
  scale_size_manual(values = point_size) +
  scale_fill_manual(values = point_colors) +
  labs(size = "Number of\nclutches observed",
       fill = "Number of\nmates observed") +
  annotate(geom = "text", y = 28, x = 16,
           label = "For visual purposes,\nonly individuals\nwith more than 2 years of\nobservations are shown",
           color = "black", size = 3, fontface = 'italic', hjust = 0)

# create the recruit plot
Rec_plot <- 
  ggplot2::ggplot() +
  geom_linerange(data = dplyr::filter(max_min_obs, age_first_cap == "J"), 
                 orientation = "y", aes(y = ring_ordered, x = min_age,
                                        xmin = min_age, xmax = max_age), 
                 size = 0.85, color = "grey60") +
  geom_point(data = dplyr::filter(eggdf_clutches, age_first_cap == "J"), 
             aes(x = est_age, y = ring_ordered, size = as.factor(n_clutches), 
                 fill = as.factor(polyandry)), shape = 21) +
  
  # geom_point(data = dplyr::filter(matings_without_egg_obs_plot, age_first_cap == "J"),
  #            aes(x = est_age, y = ring_ordered), size = 2, shape = 4) +
  scale_x_continuous(limits = c(0, 24),
                     breaks = seq(0, 24, by = 2)) +
  ylab("Locally recruited") +
  xlab("Age") +
  luke_theme +
  theme(legend.position = "none",
        axis.title.x = element_blank(),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank()) +
  scale_size_manual(values = point_size) +
  scale_fill_manual(values = point_colors)

# draw the two panels together for Fig 1
plot_of_sample_population <- 
  ggarrange(Rec_plot, Imm_plot, nrow = 2, align = "v",
            # heights = c(0.5, 0.5))
            heights = c(0.19, 0.81))

plot_of_sample_population2 <-
  (Rec_plot / Imm_plot) + 
  plot_annotation(tag_levels = 'A') + 
  plot_layout(heights = c(0.19, 0.81))

plot_of_sample_population

ggsave(plot = plot_of_sample_population2,
       filename = "products/figures/jpg/Figure_2.jpg",
       width = 6,
       height = 10, units = "in")

ceuta_egg_chick_female_data <- 
  readRDS("data/Ceuta_egg_chick_female_data.rds")

#### Modeling egg volume ----
# Full model predicting egg volume:
# Fixed effects
# - quadratic age effect (i.e., within-group deviation in age sensu van de Pol and Verhulst 2006)
# - selective appearance and disappearence (i.e., age at first and last observation)
# - female structural size (tarsus)
# - quadratic seasonal effect of first lay date (i.e., between individual effect)
# - linear within seasonal effect (i.e., within-group deviation in laydate)
# Random intercepts: nest, individual, year

# use van der Pol method to obtain the within-group deviation scores for age and
# lay date. The "x_deviation" effect is interpreted as the within individual
# effect, whereas the "first_x" (or "last_x") effect is interpreted as the between 
# individual effect.
mod_eggv_poly <-
  lmer(volume_cm ~ poly(est_age_t_deviation, 2) +
         first_age_t + last_age_t + avg_ad_tarsi + 
         laydate_deviation +
         poly(first_laydate, 2) +
         (1 | ID) + (1 | ring) + (1 | year),
       data = ceuta_egg_chick_female_data)

mod_eggv_I <-
  lmer(volume_cm ~ est_age_t_deviation + I(est_age_t_deviation^2) +
         first_age_t + last_age_t + avg_ad_tarsi + 
         laydate_deviation +
         first_laydate + #I(first_laydate^2) +
         (1 | ID) + (1 | ring) + (1 | year),
       data = ceuta_egg_chick_female_data)

# run tidy bootstrap to obtain model diagnostics
tidy_mod_eggv <-
  tidy(mod_eggv_I, conf.int = TRUE, conf.method = "boot", nsim = 1000)

model_parameters(mod_eggv_I)
random_parameters(mod_eggv_I)

# run rptR to obtain repeatabilities of random effects
rpt_mod_eggv <-
  rpt(volume_cm ~ poly(est_age_t_deviation, 2) + first_age_t + last_age_t + 
        avg_ad_tarsi + 
        laydate_deviation + poly(first_laydate, 2) +
        (1|ID) + (1|ring) + (1|year),
      grname = c("ID", "ring", "year", "Fixed"),
      data = ceuta_egg_chick_female_data,
      datatype = "Gaussian",
      nboot = 1000, npermut = 1000, ratio = TRUE,
      adjusted = TRUE, ncores = 4, parallel = TRUE)

# run partR2 on each model to obtain marginal R2, parameter estimates, and beta
# weights
R2m_mod_eggv <-
  partR2(mod_eggv_poly,
         partvars = c("poly(est_age_t_deviation, 2)",
                      "first_age_t",
                      "last_age_t",
                      "poly(first_laydate, 2)",
                      "laydate_deviation",
                      "avg_ad_tarsi"),
         R2_type = "marginal",
         nboot = 1000,
         CI = 0.95,
         max_level = 1)

R2c_mod_eggv <-
  partR2(mod_eggv_poly,
         partvars = c("poly(est_age_t_deviation, 2)",
                      "first_age_t",
                      "last_age_t",
                      "poly(first_laydate, 2)",
                      "laydate_deviation",
                      "avg_ad_tarsi"),
         R2_type = "conditional",
         nboot = 1000,
         CI = 0.95,
         max_level = 1)

# save model, tidy, rptR, and partR2 output as a list
stats_eggv_mod <-
  list(mod_I = mod_eggv_I,
       mod_poly = mod_eggv_poly,
       tidy = tidy_mod_eggv,
       rptR = rpt_mod_eggv,
       partR2m = R2m_mod_eggv,
       partR2c = R2c_mod_eggv)

save(stats_eggv_mod,
     file = "output/Stats_eggv_mod.rds")

load("output/stats_eggv_mod.rds")

model_parameters(stats_eggv_mod$mod_I, standardize = "refit")
random_parameters(stats_eggv_mod$mod_poly)
plot(allEffects(stats_eggv_mod$mod_I))

#### Libraries and data ----
source("R/project_libraries.R")
source("R/project_functions.R")
source("R/project_plotting.R")

#### Results and data----
load("output/stats_eggv_mod.rds")
ceuta_egg_chick_female_data <- 
  readRDS("data/Ceuta_egg_chick_female_data.rds")

#### Table of effect sizes (van de Pol method) ----
# Retrieve sample sizes
sample_sizes <-
  ceuta_egg_chick_female_data %>% 
  ungroup() %>% 
  summarise(Year = n_distinct(year),
            Individual = n_distinct(ring),
            Nests = n_distinct(ID),
            Observations = nrow(.))

sample_sizes <- 
  as.data.frame(t(as.data.frame(sample_sizes))) %>%
  rownames_to_column("term") %>% 
  rename(estimate = V1) %>% 
  mutate(stat = "n")

# clean model component names
mod_comp_names <- 
  data.frame(comp_name = c("Within ind. linear age",
                           "Within ind. quadratic age",
                           "Between ind. first breeding age",
                           "Between ind. last breeding age",
                           "Mother tarsus length",
                           "Within ind. lay date",
                           "Between ind. linear lay date",
                           "Between ind. quadratic lay date",
                           "Total Marginal \U1D479\U00B2",
                           "Senescence",
                           "Selective appearance",
                           "Selective disappearance", 
                           "Between ind. seasonality",
                           "Within ind. seasonality",
                           "Mother tarsus length",
                           "Total Conditional \U1D479\U00B2",
                           "Nest / Individual",
                           "Individual",
                           "Year",
                           "Residual",
                           "Nest / Individual",
                           "Individual",
                           "Year",
                           "Residual",
                           "Years",
                           "Individuals",
                           "Nests",
                           "Observations (i.e., Eggs)"))

# Fixed effect sizes (non-standardized)
fixefTable <- 
  stats_eggv_mod$tidy %>% 
  dplyr::filter(effect == "fixed") %>% 
  dplyr::select(term, estimate, conf.low, conf.high) %>% 
  as.data.frame() %>% 
  mutate(stat = "fixed")

# Fixed effect sizes (standardized)
fixef_bw_Table <- 
  stats_eggv_mod$partR2m$BW %>% 
  as.data.frame() %>% 
  mutate(stat = "fixed_bw") %>% 
  rename(conf.low = CI_lower,
         conf.high = CI_upper)

# Semi-partial R2 estimates
R2Table <- 
  bind_rows(stats_eggv_mod$partR2m$R2,
            stats_eggv_mod$partR2c$R2[1,]) %>% 
  dplyr::select(term, estimate, CI_lower, CI_upper) %>% 
  as.data.frame() %>% 
  mutate(stat = "partR2") %>% 
  rename(conf.low = CI_lower,
         conf.high = CI_upper)

# Random effects variances
ranefTable <- 
  stats_eggv_mod$tidy %>% 
  dplyr::filter(effect == "ran_pars") %>% 
  dplyr::select(group, estimate, conf.low, conf.high) %>% 
  as.data.frame() %>% 
  mutate(stat = "rand") %>% 
  rename(term = group) %>% 
  mutate(estimate = estimate^2,
         conf.high = conf.high^2,
         conf.low = conf.low^2)

# Adjusted repeatabilities
coefRptTable <- 
  stats_eggv_mod$rptR$R_boot %>% 
  dplyr::select(-Fixed) %>% 
  mutate(residual = 1 - rowSums(.)) %>% 
  apply(., 2, 
        function(x) c(mean (x), quantile (x, prob = c(0.025, 0.975)))) %>% 
  t() %>% 
  as.data.frame() %>% 
  rownames_to_column("term") %>% 
  rename(estimate = V1,
         conf.low = `2.5%`,
         conf.high = `97.5%`) %>% 
  mutate(stat = "RptR")

# Store all parameters into a single table and clean it up
allCoefs_mod <- 
  bind_rows(fixef_bw_Table,
            R2Table,
            ranefTable, 
            coefRptTable, 
            sample_sizes) %>% 
  bind_cols(.,
            mod_comp_names) %>%
  mutate(coefString = ifelse(!is.na(conf.low),
                             paste0("[", 
                                    round(conf.low, 2), ", ", 
                                    round(conf.high, 2), "]"),
                             NA),
         effect = c(rep("Fixed effects \U1D6FD (standardized)", nrow(fixef_bw_Table)),
                    rep("Partitioned \U1D479\U00B2", nrow(R2Table)),
                    rep("Random effects \U1D70E\U00B2", nrow(ranefTable)),
                    rep("Adjusted repeatability \U1D45F", nrow(coefRptTable)),
                    rep("Sample sizes \U1D45B", nrow(sample_sizes)))) %>%
  dplyr::select(effect, everything())

# re-organize model components for table
allCoefs_mod <-
  allCoefs_mod[c(5, 1:4, 6:9, 16, 15, 10:14, 17:28), ]

# draw gt table
eggv_mod_table <- 
  allCoefs_mod %>% 
  dplyr::select(effect, comp_name, estimate, coefString) %>% 
  gt(rowname_col = "row",
     groupname_col = "effect") %>% 
  cols_label(comp_name = html("<i>Egg volume</i>"),
             estimate = "Mean estimate",
             coefString = "95% confidence interval") %>% 
  fmt_number(columns = vars(estimate),
             rows = 1:24,
             decimals = 2,
             use_seps = FALSE) %>% 
  fmt_number(columns = vars(estimate),
             rows = 25:28,
             decimals = 0,
             use_seps = FALSE) %>% 
  fmt_missing(columns = 1:4,
              missing_text = "") %>% 
  cols_align(align = "left",
             columns = vars(comp_name)) %>% 
  tab_options(row_group.font.weight = "bold",
              row_group.background.color = brewer.pal(9,"Greys")[3],
              table.font.size = 12,
              data_row.padding = 3,
              row_group.padding = 4,
              summary_row.padding = 2,
              column_labels.font.size = 14,
              row_group.font.size = 12,
              table.width = pct(60))

eggv_mod_table